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INTRODUCTION 


Modern  digital  ionoaondes  provide  many  new  opportunities  in  ionospheric 
observations  and  promise  significant  progress  in  understanding  the  dynamic 
processes  in  the  ionosphere.  The  capabilities  and  sensitivity  of  digital 
ionoaondes  will  be  demonstrated  by  using  the  NOAA  ionosonde  (Grubb,  R.  N., 
1979)  as  a  prototype  example.  The  emphasis  will  be  on  the  mathematical 
principles  of  ionogram  processing,  which  will  be  applicable  to  most  digital 
ionograms  depending  on  their  design  and  mode  of  operation.  At  the  same  time, 
this  report  may  give  some  guidance  for  the  design  of  new  digital  ionoaondes, 
depending  on  the  user's  needs. 


THE  DATA 


Ionograms  obtained  with  classical  ionoaondes  produced  a  display  of  the 
travel  times  or  virtual  heights  of  echoes  as  a  function  of  frequency.  Both 
the  frequency  and  the  virtual  heights  were  defined  in  some  analog  fashion  and 
were  of  limited  accuracy.  Modem  digital  ionoaondes  derive  the  frequencies 
from  a  frequency  synthesizer  with  very  high  precision  and  also  measure  the 
amplitude  and  phase  of  each  echo  in  addition  to  the  virtual  heights.  As  we 
will  see  later,  it  is  especially  this  phase  which  provides  new  and  more 
accurate  information.  In  a  typical  mode  of  operation,  the  NOAA  ionosonde  has 
four  short  dipoles  as  receiving  antennas  located  in  the  comers  of  a  square 
and  oriented  along  the  diagonals  of  the  square.  A  two-channel  quadrature 
receiver  is  connected  to  one  or  other  diagonal  antenna  pair  during  one 
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transmit- receive  cycle.  A  full  data  set  is  obtained  by  tranoaitting  four 
pulses,  two  at  a  frequency  f  and  two  at  a  slightly  higher  frequency  f  +  A  f, 
where  Af  is  kept  constant  through  the  frequency  sweep  and  is  usually  8  kHz. 
Such  a  set  then  gives  virtual  heights,  the  quadrature  components  equivalent  to 
amplitude  and  phase  for  each  echo  for  four  antennas  and  two  frequencies. 
Since  the  echo  detection  scheme  is  based  on  lo-ps  samples,  and  differences 
in  travel  time  of  the  eight  echoes  in  a  set  are  email  compared  with  the 
sampling  interval,  care  is  taken  that  the  echo  information  is  recorded  at  the 
same  reference  time  relative  to  the  pulse  transmission  closest  to  the  average 
maximum  amplitude. 

Phase  measurements  are  accurate  to  approximately  1°  if  the  signal-to- 
noise  ratio  is  high.  The  quadrature  components,  the  output  of  the  linear 
receiver  channels,  are  digitized  to  12-bit  numbers  (including  sign).  With 
each  transmitted  pulse,  a  delayed  calibration  pulse  is  injected  into  the 
preamplifiers  at  the  center  of  the  antenna  square,  is  recorded,  and  permits 
correction  of  possible  differences  in  cable  length  from  antenna  to  receiver, 
differences  in  the  preamplifier  characteristics,  and  differences  between 
receiver  channels.  The  noise  level  is  also  recorded,  averaged  over  a  window 
between  the  transmitted  pulse  and  the  lowest  expected  E-region  echo.  This 
way,  echoes  with  poor  signal- to-noise  ratio  can  be  eliminated  in  the  data 
analysis. 

THE  INFORMATION  CONTENT  OF  THE  PHASE 

The  phase  value  $  obtained  for  an  individual  echo  is  linearly  related  to 
the  phase  height  (or  phase  integral) .  The  constant  term  in  this  relation  is 
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unknown  and  normally  cannot  be  determined.  However,  the  relation  between  the 
variation  of  the  phase  with  time,  frequency,  location,  and  Ionospheric 
parameters  Is  well  known.  This  can  be  expressed  in  the  following  way: 


34  3  A  34  34 

d*  “  3x  ^  +  3y  dy  +  3t  dt  +  3f  df  +  d*o*  ( 1} 


$ 

■i 
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(It  is  assumed  that  the  four  receiving  antennas  are  in  one  horizontal  plane: 

otherwise,  another  term  for  the  height  variation  has  to  be  included.)  In  this 

equation,  the  spatial  variation  of  the  phase  yields  the  angle  of  arrival  of 

the  echo;  the  temporal  variation  yields  the  Doppler  frequency  or  Doppler 

velocity  (change  of  phase  height  with  time):  and  the  variation  with  frequency 

gives  a  more  accurate  estimate  of  the  virtual  height  according  to  the 

principle  of  stationary  phase.  The  quantity  d4  accounts  for  the  difference 

o 

in  antenna  orientation. 

We  assume  now  that  one  antenna  pair  is  in  the  north-south  direction  and 
the  other  in  the  east-west  direction,  and  label  the  antennas  N,  E,  S,  W.  We 
further  use  a  cartesian  coordinate  system  with  the  x-axis  pointing  east  and 
the  y-axis  pointing  north,  and  its  center  in  the  center  of  the  antenna  square. 
To  separate  the  several  phase  variations,  we  use  the  following  antenna 
frequency  switching  sequence: 
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(More  sophisticated  antenna  frequency  switching  sequences  were  suggested  by 
Wright,  Pitteway  et  al,  1982).  We  now  assume  that  D  is  the  length  of  the 
diagonal  of  the  antenna  square  and  At  the  time  interval  between  consecutive 
pulses.  It  is  now  easy  to  derive  the  following  set  of  equations: 
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Each  of  the  unknowns  in  this  set  of  equations  can  be  obtained  by  two  indepen¬ 
dent  combinations  of  two  or  four  measured  phase  values.  The  simplest  solution 
is  bo  take  the  mean  value  as  the  result  for  each  unknown  and  the  mean 
difference  as  an  error  estimate.  The  same  mean  values  are  obtained  in  a 
least-squares  fit  if  we  consider  the  measured  phase  values  as  deviations  from 
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a  fictitious  phase  value  for  the  center  of  the  4-dimensional  box  (x,  y,  t,  f ) . 


Before  discussing  the  phase  deviations  in  more  detail,  it  is  important  to 


remember  the  relation  between  the  phase  $  and  the  phase  path  p*.  The  phase 


path  p*  is  defined  by: 


where  y  is  the  phase  refractive  index  and  ds  the  path  element.  On  the  other 


hand,  we  have: 


P*  -  «♦ 


.  _  2ir 
♦  “  e  f 


ANGLE  OF  ARRIVAL 


Outside  the  ionosphere,  y=l  and  grad  p*  according  to  (3)  is  a  unit  vector 


pointing  in  the  direction  of  propagation;  in  other  words,  the  z-component  of 


grad  p*  is  negative  for  an  echo  reflected  from  the  ionosphere.  We 


obtain  and  — from  (2)  and  (4)  and  have,  therefore, 
3x  3y 


hi.  .  .  M  ,\W 

dz  1  '^X-'  'W  ’  1  * 


The  direction  toward  the  apparent  reflection  point  in  the  ionosphere  is  then 


given  by  a  unit  vector  £  with 


v  ‘  (V  V  “  “  grad 


S 


The  zenith  distance  £  and  the  azimuth  angle  u  (measured  from  east  over  north 


consistent  with  the  cartesian  coordinate  system)  are  then  given  by 


DOPPIER  FREQUENCY  AND  VELOCITY 


Front  (2),  we  obtain  the  change  of  phase  with  time  and  can  define  a 
Doppler  frequency  by 


*  1  3£ 

fa  =  n  a- 


(10) 


In  many  cases,  it  may  be  more  interesting  to  define  a  Doppler  velocity  by 

V  =|£l=X_|i  (11) 

d  3t  2u  3t 

It  must  be  pointed  out  that  in  the  above  definition  is  the  total  path 
change  with  time.  In  case  of  echoes  reflected  from  a  layer  moving  as  a  whole, 
one-half  the  velocity  defined  above  would  be  the  true  velocity  of  the  motion 
of  the  layer.  In  this  case,  the  velocity  would  be  constant  for  all  frequen¬ 
cies  reflected  from  such  a  moving  layer,  while  the  Doppler  frequency  would 
increase  proportional  with  the  radio  frequency.  Since  the  phase  path  is 


increasing  when  the  distance  to  the  layer  increases,  the  Doppler  velocity  will 
be  positive  and,  correspondingly,  the  velocity  will  be  negative  when  the  dis¬ 
tance  decreases.  The  interpretation  of  such  a  velocity  is  not  always  simple, 
flor  example,  assume  that  the  electron  density  locally  increases  somewhere 
along  the  path  without  any  change  near  or  at  the  reflection  level.  In  this 
case  the  phase  path  decreases,  the  group  path  increases,  and  the  true  path 
length  is  unchanged  (Paul  et  al,  1974). 

A  typical  example  of  this  type  of  measurement  is  shown  in  figure  1.  The 
data  were  taken  during  the  sunrise  period  on  8  January  1981  at  Brighton, 
Colorado,  and  the  results  shown  were  derived  as  averages  over  three  consecu¬ 
tive  frequencies  near  the  F-region  peak.  Ordinary  and  extraordinary  echoes 
were  analyzed  separately  at  corresponding  frequencies  and  are  marked  by  o  and 
x  in  the  plot.  The  echoes  came  approximately  from  southeast  as  expected,  and 
approximately  10°  off  vertical.  The  Doppler  velocity  of  -50  m/s  is  of  course 
related  to  the  velocity  with  which  the  ionization  front  approaches.  The 
lowest  part  of  figure  1  shows  the  deviation  of  the  MUF  (3000)  from  the  rapid 
linear  increase  of  this  quantity  (to  be  discussed  later),  indicating  some 
disturbance  by  gravity  waves  which  is  also  visible  in  the  temporal  variation 
of  the  other  quantities  shown.  The  close  agreement  between  ordinary  and 
extraordinary  results  indicates  a  high  reliability  of  the  data  and  their 
analysis. 

GROUP  PATH 

According  to  the  principle  of  stationary  phase,  the  group  path  p'  and 
phase  path  p*  are  related  in  the  following  ways 


df  -  2  it  df 


(12) 


where  c  is  the  velocity  of  light. 

For  vertical  reflection,  the  corresponding  relation  between  virtual 
height  h  and  phase  height  h*  is  valid.  With  (12)  we  have  the  possibility  of 
measuring  the  virtual  height  or  group  path  independently  from  the  measurement 
of  the  travel  time.  Obtaining  the  travel  time  means  first  selecting  a  certain 
characteristic  of  the  shape  of  the  transmitted  pulse  (leading  edge,  steepest 
slope,  maximum  amplitude,  etc)  which  is  also  expected  to  be  found  in  the 
reflected  echo.  The  time  difference  between  the  appearance  of  this  character¬ 
istic  in  the  transmitted  and  reflected  pulse  is  then  the  travel  time.  This 
procedure  has  limited  accuracy,  since  the  pulse  shape  can  be  distorted  between 
transmission  and  reception  (Paul,  1979).  in  general,  the  phase  is  less  sensi¬ 
tive  to  distortion  and  the  phase  change  over  a  small  frequency  interval  gives 
a  more  accurate  estimate  of  the  group  path  or  virtual  height.  The  phase  may 
change  by  several  full  rotations  plus  a  fraction  of  a  rotation  over  a  given 
frequency  interval.  Only  the  fraction  is  obtained  by  the  phase  comparison, 
and  the  number  of  full  rotations  has  to  be  determined  from  the  travel  time. 

The  frequency  increment  Af  =  8  kHz  corresponds  to  a  basic  virtual  height 
unit  of  18.75  km  for  one  full  phase  rotation.  If  phase  changes  over  this 
frequency  interval  can  be  measured  with  an  accuracy  of  1°,  the  resulting 
virtual  heights  have  an  accuracy  of  about  50  ra.  For  a  smaller  Af  and  the  same 
phase  resolution,  the  accuracy  of  the  virtual  heights  would  decrease.  On  the 
other  hand,  for  a  virtual  height  of  300  km,  the  phase  changes  by  16  full  rota- 
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tions  over  the  frequency  interval  of  8  kHz  and  a  much  larger  value  of  Af  may 
cause  difficulties  in  determining  the  nutber  of  full  rotations  unanbiguously. 
In  addition,  for  a  large  Af  the  linear  approximation  of  the  phase  variation  in 
the  frequency  becomes  inaccurate  and  a  higher-order  function  for  $  (f)  would 
have  to  be  used  in  the  analysis. 

POLARIZATION  DISCRIMINATION 

Assume  that  ordinary  and  extraordinary  echoes  at  a  given  frequency  are 
well  separated  in  virtual  heights.  The  phases  of  the  north-south  antenna  pair 
then  should  be  different  from  phases  at  the  east-west  pair  by  d$^  =  +  90°  or 
d<£  =  -  90* ,  depending  on  the  sense  of  rotation  [see  equation  (1)].  For  the 

purposes  of  identifying  an  echo  as  ordinary  or  extraordinary  it  is  not 
necessary  to  solve  for  d$^  .  Comparing  the  mean  phase  values  of  the  north- 
south  pair  with  the  mean  phase  value  of  the  east-west  pair,  we  have  the  one  or 
other  component,  depending  on  whether  the  phase  difference  is  close  to  +90°  or 
close  to  -90°.  This  is  still  true  if  the  pulse  is  reflected  obliquely. 

SEPARATION  OF  OVERLAPPING  ORDINARY  AND  EXTRAORDINARY  ECHOES 

In  most  ionograms,  there  are  some  frequency  ranges  where  ordinary  and 
extraordinary  echoes  have  approximately  equal  travel  times  and  the  tvo  traces 
cross  each  other.  If  the  two  echoes  overlap  sufficiently,  we  obtain  an  inter¬ 
ference  pattern  as  a  function  of  frequency  and  the  distinction  between  the  tv» 
modes  in  terms  of  travel  time  is  no  longer  possible.  The  resulting  sequence 
of  maxima  and  minima  as  a  function  of  frequency  is  not  in  phase  at  the  four 
receiving  antennas  except  that  both  echoes  are  returned  exactly  vertically. 
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If  amplitude  and  phase  at  each  of  the  four  antennas  are  recorded  for  a  given 
frequency,  the  separation  of  the  two  echoes  is  possible  under  rather  general 
conditions. 

It  is  assumed  that  the  two  components  are  circularly  polarized  with 
opposite  sense  of  rotation.  They  will  appear  to  be  elliptically  polarized,  if 
they  are  not  reflected  vertically.  As  mentioned  earlier,  the  apparent  direc¬ 
tion  to  the  reflection  point  can  be  described  by  the  zenith  distance  £  and  the 
azimuth  angle  a  and  the  components  of  the  unit  vector  w  in  (7)  are  then  given 
by: 


w  =  [sin  £  cos  a,  sin  £  sin  a,  cos  £] . 

The  electric  vector  is  perpendicular  to  w,  and  a  unit  vector  e  in  the 
direction  of  the  electric  vector  must  satisfy  the  following  conditions: 

e  +  e  +  e  =1  '14) 

x  y  z 

and  ew+ew+ews0.  (15) 

xx  y  y  z  z 

At  an  instant  when  the  electric  vector  is  in  the  (y,  z)  plane,  we  have  with 
(13),  (14),  and  (15), 

e  »  [0,  cos  C/sc,  sin  £  sin  o/sc]  (15) 

N 


with 


sc  *  (1  -  sin^  £  cos^  a)^. 

When  e  is  in  the  (x  ,  z)  plane,  we  have 

eE  *  *cos  ^8S»  0  •  8in  ^  cos  a/ss] 


(17) 


(18) 


with 


2  2  Vi 

ss  =*  (1  -  sin  £  sin  a)  . 


(19) 


The  amplitude  observed  at  the  north  and  south  antenna  (oriented  in  the  north- 
south  direction)  would  then  be  proportional  to  e^  in  (16)  and  the  amplitude 
observed  at  the  east  and  west  antenna  would  be  proportional  to  e^  in  (18)  with 
the  same  proportionality  factor.  In  other  words,  the  amplitude  ratio  r 
between  the  north-south  and  the  east-west  antenna  pairs  is  given  by: 


r  *  ss/sc. 

A 


(20) 


Corresponding  to  the  off  vertical  incidence  of  an  echo,  we  observe  a  phase 

difference  between  the  north  and  south  antenna  and  also  between  the  east  and 

west  antenna.  Those  phase  differences  are  proportional  to  e^  in  (16)  and  (18) 

respectively.  Defining  A$  as  the  phase  difference  between  the  north  antenna 

N 

and  a  ficticious  antenna  in  the  center  of  the  array,  and  correspondingly  a 
phase  difference  AA  ,  we  have 

u 


A$ 


sin  £  sin  a/sc 


(21) 


and 


A* 


-IT 


sin  £  cos  a/ 


ss 


(22) 
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where  X  is  the  radio  wavelength  and  D  is  the  diagonal  distance  between  antenna 
pairs  as  in  (2) . 

We  now  can  write  the  equations  for  the  quadrature  components  of  two  over¬ 
lapping  echoes  with  their  electric  vectors  rotating  in  opposite  directions, 
neglecting  temporal  variation  over  the  interval  of  two  consecutive  pulses*  We 

use  the  symbols  a,  e  ,  b,  c  for  amplitude  and  phase  difference  on  the  north 

N  E 

and  east  antenna  for  one  component  and  b,  6  ,  d,  6  correspondingly  for  the 

other  component.  The  phases  of  the  two  components  for  the  center  of  the  array 

are  described  by  <b  and  $  .  The  indices  X  and  Y  Indicate  the  cosine  and  sine 

1  2 

term  of  the  quadrature  component,  while  the  symbols  N,  E,  S,  and  W  are  used 
for  the  magnitude  of  the  quadrature  components  at  the  four  antennas. 


W„ 


=  a  cos  (♦1  +  e^)  +  b  cos  ($2  +  <5^) 

*  a  sin  <  4>  +  e  )  +  b  sin  ($  +  6  ) 

IN  2  N 

=  -c  sin  +  e  )  +  d  sin  ($  +  6  > 

IE  2  E 

=  c  cos  (<fr  +  e  )  -  d  cos  ($  +  6  ) 

1  E  2  G 

=  a  cos  ($  "  £  )  +  b  cos  ($  -  6  ) 

IN  2  N 

*  a  sin  (♦,  -  e  )  +  b  sin  ($  -  5  ) 

IN  2  N 

*  -c  sin  -  e  )  +  d  sin  -  6  ) 

IE  2  E 


(23) 
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W  -  c  cos  (A  -  e  )  -  d  cos  (A  -  i  ) . 
y  1  E  2  E 

These  ere  eight  equations  for  10  unknowns.  In  addition,  we  have  according 
to  (20), 

—  *  ss./sc,  (24) 

C  II 

where  ss^  and  sc^  are  the  quantities  defined  in  (17)  and  (19)  for  €^and  , 
the  zenith  distance  and  azimuth  for  one  component;  according  to  (21)  and  (22), 

CN  "  “*  X  Si"  S  Sin  a^8C\ 

(25) 

D 

e  =  -7t  —  sin  E  cos  a  /s s„ . 

E  A  1  11 

The  corresponding  set  of  equations  for  the  other  component  is  then 

I  "  SS2/8C2  (26) 

5N  “  I  sin  C2  sin  <*2/sc2 

(27) 

$E  *  J  ®in  C2  cos  a2/ss2 . 

Equations  (24)-(27)  add  six  more  equations  for  four  more  unknowns  and  we  have 
now  14  equations  for  14  unknowns.  Unfortunately,  those  equations  are  non¬ 
linear  and  do  not  offer  an  easy  exact  solution,  if  such  an  exact  solution 
exists  at  all.  In  the  following,  we  will  show  how  an  approximate  solution  can 
be  found  which  can  be  improved  by  iteration. 
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First,  we  define  eight  combinations  of  the  quadrature  components 


R  =  N  +  S  =  2  a  cos  cos  £  +  2  b  cos  ♦  coe  5 

1  X  X  IN  2  N 

R  =  N  -  S  =  -2  a  sin  $ .  sin  e  -  2  b  sin  $  sin  5 

2  X  X  I  N  2  N 


R->  =  +  =  2  a  sin  cos  £  +  2  b  sin  $  cos  4 

3  Y  Y  1  N  2  N 

R4  =  =  2  a  cos  ^  sin  +  2  b  cos  ^  sin  6^  (28) 


R„  =  E  +  W  =  -2  c  sin  4.  cos  C  +  2  d  sin  <k  cos  6 
5  X  X  1  E  2 


R,  =  E„  -  w„  =  -2  c  cos  sin  e  +  2  d  cos  $  sin  6 
O  X  X  ■  12s  2  IS 


R7  =  =  2  c  cos  4^  cos  ££  -  2  d  cos 


cos 


6 

E 


E  -  W  =  -2  c  sin  $  sin  e  +  2  d  sin  ♦  sin  6  . 
Y  Y  IE  2  E 


We  now  use  the  following  identities 


2  a  =  a  *  c  +  (a-c) 


2  c  -  a  +  c  -  (a-c) 


2  b  ■  b  +  d  +  (b-d) 
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2d*b+d-(b-d) 


and  define 


T  =  R  -  (a-c)  cos  4  cos  e  -  (b-d)  cos  ♦  cos  6 
11  IN  2  N 


T2  *  R2  +  (a-c)  sin  sin  +  (b-d)  sin  $2  sin  6^ 


T  =  R  -  (a-c)  sin  A  cos  e  -  (b-d)  sin  cos 
3  J  IN  2  N 


T  =  R  -  (a-c)  cos  4  sin  e  -  (b-d)  cos  i  sin  5 
4  4  Y1  N  2  N 


T  =  R  +  (a-c)  sin  4>  cos  e  -  (b-d)  sin  $  cos  6 

D  J  I  E  4  E 


T_  *  R  +  (a-c)  cos  sin  e  -  (b-d)  cos  sin  S 
ob  IE  2  E 


T_  =  R  -  (a-c)  cos  4  cos  e  +  (b-d)  cos  $  cos  6 
7  7  1  E  2  E 


T  =  R  +  (a-c)  sin  4  sin  e  -  (b-d)  sin  $  sin  6  • 
8  8  1  E  2  E 


According  to  (28),  the  equation  system  (29)  is  equivalent  to 


T  *  (a+c)  cos  4  cos  e  +  (b+d)  cos  4  cos  6 
1  IN  2  N 


T  *  -(a+c)  sin  4  sin  e  -  (b+d)  sin  4  sin  6 
2  IN  2  l 


T  *  (a+c)  sin  $  cos  e  +  (b+d)  sin  cos  6 
3  I  N  2  N 


(29) 
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(30) 


*  (a+c)  cos  sin  +  (b+d)  cos  $2  8^n 
Tg  =*  -(a+c)  sin  cos  +  (b+d)  sin  $2  cos  $E 

=*  -(a+c)  cos  sin  +  (b+d)  cos  8^n  ^g 

=  (a+c)  cos  cos  £g  —  (b+d)  cos  $2  008  ^g 

T  =  -(a+c)  sin  sin  e_  +  (b+d)  sin  $  sin  6  . 

8  i  ® 

Using  t'ne  symbols 


S1  =  sin  4»1 
=  cos  4>1 

V  sin  <)»2 

C  =*  COS  4> 

2  2 

we  derive  from  (30) 

S  T  -C  T  =  -(b+d)  sin  ($  -  $  )  cos 
1113  *  1  ” 

S  T  +  C.  T_  *  -(b+d)  sin  (<fr  “  ♦-)  8^n 
14  12  z  1 

S1  T7+  C1  T5  *  (b+d>  8in  (*2  ”  *1*  008  SE 
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(31) 


•  Yrr  ~ f  ii 


«• 


S.  “  C,  T  -  -(b+d)  sin  (<j>  -  4  >  sin  6 

1  o  l  o  <  1  E 

and 


S2  T1  “  C2  T3  *  (a+c>  8in  ^ 2  ~  008  EN 

S,  T  +  C,  T  -  (a+c)  sin  (<j>  -  ♦  )  sin  e  (32) 

2  4  2  2  2  1  N 

S2  T?  +  C2  T5  =  <a+c)  8in  c°s  eE 

S2  T&  -  C2  Tg  =  -(a+c)  sin  <$2  -  sin  e^.. 

We  now  obtain  from  (31) 


(S1  T1  -  C1T3)2  +  (S1  T4  +  C1  T2)2  =  (S1T?  +  Ct  Tg)2  +  (S^  -  Tg)2  03) 

and  a  similar  equation  from  (32)  with  the  only  difference  that  and  are 
replaced  by  S2  and  C2  respectively.  In  general,  we  can  write  instead  of  (33) 

S2(T12+T42-T62-T72)-2SC(T1T3-T2T4+T5T7-T6T8)+C2  <T22'Pr32_T52”T82)  “  0  04) 

If  we  know  the  T  ,  i  *  1,  8  and  assume  that  C  *  0  we  can  derive  from  (34)  a 
quadratic  equation  for  S/C  and  its  roots  are  tg  and  tg  $2« 

We  see  now  how  an  iteration  process  can  be  formed.  In  a  first  step,  we 
assume  that  the  zenith  angles  are  relatively  anall 


sin2  <  1  i  -  1,  2 


05) 
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u 


which  means,  according  to  (17),  (19),  and  (20),  that 

a  ID  c  and  b  X  d  ; 

in  other  words  we  neglect  terms  with  factors  (a  -  c)  and  (b  -  d)  in  (29),  and 
have 


i  —  1,  2,  • • • • 8 , 

and  obtain  tg  and  tg  $2  from  (34)  Unfortunately  this  means  an  ambiguity 
in  $  and  by  ±  it.  This  will  be  discussed  later.  Knowing  and  we 
obtain  6  from  the  first  pair  of  equations  in  (31),  6  from  the  second  pair, 
and  b  +  d  from  either  one  pair.  Correspondingly  we  obtain  e^,  and  a  +  c 
from  (32). 

Now  assume  that  we  change  <J>  1  by  it.  We  see  that  this  would  cause  a  change 
of  the  sign  on  the  left  and  right  side  of  all  equations  in  (31),  and  there¬ 
fore  6  and  6  remain  unchanged;  while  in  all  equations  (32),  only  the  right 
N  E 

side  would  change  the  sign  which  has  to  be  compensated  by 

changing  e  and  e  by  it.  Similarly,  a  change  of  $  by  it  would  cause  a 

change  of  6  and  6  by  if  and  leave  e  and  e  unchanged.  If  $  and  ♦  are 
N  E  N  E  12 

changed  by  it,  all  four  phase  differences  6  ,  5  ,  e  ,  and  e  have  to  be  changed 

N  E  N  E 

by  it.  With  the  reasonable  restriction  that  6  ,  6  ,  e  ,  and  e  are  all  to  be 

N  E  N  E 

within  t  it/2,  we  can  eliminate  in  most  cases  the  ambiguities  in  ^  and 
This  restriction  is  consistent  with  the  assumption  of  snail  zenith  distances, 
at  least  for  low  frequencies.  For  higher  frequencies,  D/X  >1,  the 
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quantities  6  ,  6  ,  e  ,  and  e  can  become  relatively  large  even  for  small 
NEE  G 

values  of  £;  here,  continuity  with  frequency  can  be  used  to  resolve  the 
ambiguities. 


Knowing  e  ,  e 
N  E 


u  * 


V  = 


w  = 


and  a  +  c  and  using 


X_ 

EN  UD 


e  — 
E  UD 


a  +■  c 


(36) 


we  obtain  with  (24)  and  (25), 


sin  *  (u2  +  v2  -  2u2v2),/2/  (1  -  u2v2  ft  (37) 

sin  <*t  «  -u  (1  -  v2)1/a/  (u2  +  v2  -  2u2v2)1/z  (38) 

cos  <*1  *  -v  (1  -  u2)^2/  (u2+  v2  -  2u2v2)^2  (39) 

a  =  w  (1  -u2)1/2/  (  (1  -  v2)1/2+  (1  -  u2)1/25  (40) 

c  -  w  (1  -v2)1/2/  £  (1  -  v2)1/2+  (1  -  u2)1/2l.  (41) 


Exchanging  e  by  6  ,  e  by  6  ,  a  by  b,  c  by  d,  and  the  index  1  by  2,  we  obtain 
N  N  E  E 

the  corresponding  parameters  for  the  other  component  from  (37)-(41). 
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We  now  have  a  first-order  estimate  of  all  unknowns  and  can  improve  our 


estimates  of  the  T.,  by  subtracting  the  respective  terms  from  R^  according  to 
(29),  and  the  next  iteration  can  be  computed.  This  process  was  tested  on  many 

simulated  and  real  data  and  was  found  to  be  rapidly  convergent.  Very  often, 

the  first  iteration  was  already  sufficient  to  obtain  reasonably  accurate 
results.  After  convergence  of  this  iteration  process  has  been  accomplished 

for  a  frequency  f,  we  apply  the  same  procedure  to  the  data  at  f  +  Af.  After 

the  phases  for  each  component  at  f  and  f  +  Af  are  computed,  their  virtual 
heights  are  then  obtained  from  the  phase  difference  over  the  frequency  inter¬ 
val  Af.  An  example  of  this  type  of  analysis  is  demonstrated  in  figure  2.  The 
symbol  M  (mixed)  indicates  the  virtual  heights  observed  by  measuring  the 
travel  time  of  the  overlapped  echoes,  while  the  dots  (ordinary  component)  and 
the  crosses  (extraordinary  component)  are  the  results  of  the  mode  separation 
process.  It  should  be  kept  in  mind  that  the  nominal  pulse  duration  is  60  Ps 
and  the  actual  pulse  duration  at  the  receiver  output  about  100  Ps, 
corresponding  to  a  pulse  length  on  the  virtual  height  scale  of  9  km  and  15  km, 
respectively. 


THE  MUF ( 3000) 

The  maximum  usable  frequency  over  a  3000  km  path  MUF(3000)  and  its  ratio 
to  the  ordinary  critical  frequency  foF2,  the  propagation  factor  M(3000)  is  one 
of  the  standard  parameters  in  scaling  ionograms.  Instead  of  the  table 
published  by  URSI  (Piggot,  W.R. ,  and  K.  Rawer,  1972),  we  use  an  empirical 
algebraic  expression  for  the  MUF  function  which  fits  the  tabulated  values 
closely  (better  than  0.5%  in  the  range  200  <  h  <  600  km).  We  define  the  MUF 
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function 


Normally,  the  MUF(3000)  is  obtained  with  very  high  precision,  better  than  foF2 
or  M( 3000 ) .  Besides  its  importance  for  hf  communication,  the  MDF(3000)  is 
very  sensitive  to  natural  or  artificial  disturbances  in  the  F-region  and  can 
be  used  as  an  indicator  for  such  phenomena.  In  the  numerical  process,  r  is 
evaluated  for  a  sequence  of  virtual  heights-frequency  pairs  and  the  running 
average  over  seven  values  is  computed.  When  this  average  starts  to  decrease 
with  increasing  frequency,  a  least  squares  parabolic  fit  over  the  last  seven 


individual  values  of  r  is  used  to  determine  r  ,  which  is  equal  to  the 

max 

MUF(3000).  Figure  3  shows  a  sequence  of  MUF(3000)  values  computed  from 
ionograms  taken  in  2-minute  intervals  over  a  period  of  about  8  hours.  The 
accuracy  and  sensitivity  of  this  quantity  can  better  be  seen  in  the  insert, 
where  a  portion  of  the  same  data  is  displayed  on  a  5- time 3- expanded- frequency 
scale.  (Note  that  each  point  is  obtained  independently  from  one  ionogram,  and 
no  temporal  snoothing,  averaging,  or  filtering  was  applied.) 

F  REGION  PEAK  PARAMETERS 

Perhaps  the  most  important  parameter  obtained  by  vertical  sounding  of  the 
ionosphere  is  the  penetration  frequency  or  critical  frequency  in  the  ordinary 
mode  foF2.  Actually,  it  is  obtained  by  extrapolation  when  the  virtual  heights 
as  a  function  of  frequency  approach  infinity.  As  always,  an  extrapolate* 
quantity  is  likely  to  be  less  accurate  than  a  direct  mea^ur  absent  oi  inter¬ 
polated  quantity.  This  has  to  be  kept  in  mind  when  peak  parameters,  like 
height  of  maximum  or  scale  height,  are  discussed.  Even  if  very  precise 
virtual  height  and  frequencies  are  available,  small  disturbances  near  the  F- 
region  peak  can  cause  relatively  large  errors.  It  is  generally  assuned  that 
the  F-region  peak  can  be  best  described  by  a  parabolic  model.  Under  this 
assumption,  it  has  been  shown  (Paul  et  al,  1981)  that  consecutive  frequency 
differences  divided  by  the  corresponding  virtual  height  differences  decrease 
linearly  as  a  function  of  frequency  and  their  extrapolation  to  zero  gives  the 
critical  frequency.  This  procedure  can  be  used  for  ordinary  and  extraordinary 
mode.  We  obtain  the  remaining  peak  parameters,  the  height  of  the 

maximum  h  and  the  half  thickness  y  (twice  the  scale  height)  in  the  following 
m 


way.  We  first  calculate  parabolic  coefficients  pc  for  the  last  few  frequen- 


cies  below  the  critical  frequency  where  virtual  heights  were  observed.  Those 
coefficients,  of  course,  depend  on  the  critical  frequency.  We  then  obtain 
equations  of  the  form 

h  .  *  h  +  (pc  -1)  y  (47) 

i  m  *1  m 

Applying  the  method  of  least  squares  to  several  (5-10)  of  those  equations, 

we  obtain  h  and  y  .  An  improved  estimate  of  the  three  peak  parameters  can 
m  m 

be  obtained  in  an  iterative  least- squares  process  in  which  we  linearize  the 
dependence  of  pc  on  the  critical  frequency,  which  cam  only  be  done  over  very 
small  frequency  intervals  (<50  kHz).  This  can  be  important  when  the  peak  is 
slightly  distorted  (eg,  by  gravity  waves).  If  fc  is  the  critical  frequency 
(ordinary  or  extraordinary  mode),  we  write  instead  of  (47) 

,  dpc, 

h  .  =  h  +  (pc.  +  ■  Afc  -  1)  y  (48) 

i  m  dfc  •'m 

for  the  three  unknowns  hm,  ym,  and  Afc.  The  quantity  Afc  is  now  a  correction 
to  the  critical  frequency,  and  the  process  can  be  repeated  with  the  new 
critical  frequency.  The  algorithm  for  computing  the  parabolic  coefficients 
will  be  outlined  later. 

This  method  has  been  applied  to  hundreds  of  ionograms  and  also  was  tested 
with  simulated  data.  In  most  cases,  very  good  agreement  between  these  results 
and  full  profile  computations  was  found  except  in  situations  where  systematic 
differences  were  expected,  mainly  when  a  very  pronounced  FI  layer  was  present 
and  its  critical  frequency  was  relatively  close  to  foF2.  In  those  cases, 
empirical  corrections  depending  on  the  ratio  foFl/foF2  can  be  applied. 
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t- 


Results  of  those  studies  trill  be  reported  later.  If  this  peak  fitting  process 
is  applied  to  ordinary  and  extraordinary  components  separately,  the 
differences  of  the  results  give  some  indication  of  the  reliability  of  the 
computed  parameters. 

The  peak  fitting  program  together  «rith  the  MCJF(3000)  computation  can  be 
considered  as  a  first  order  automatic  ionogram  processing  routine  and  has  been 
applied  directly  to  many  tapes  with  sequences  of  500  and  more  ionograms 
without  interruptions. 


The  parabolic  coefficients  can  be  computed  in  the  following  way.  The 
virtual  height  h  can  be  expressed  as 


h 


(49) 


where  y  is  the  group  refractive  index  and  h^  is  the  reflection  height.  The 
principle  is  now  to  write  the  identity 


h 


dh 


(50) 


where 


fN2/f 


2 


and  fN  is  the  plasma  frequency. 
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The  product  p  /  1-X  Is  the  ratio  of  the  field  included  group  refractive 
index  to  the  no-magnetic-field  group  refractive  index,  which  is  easy  to  eval¬ 
uate  and  usually  is  close  to  unity.  If  a  closed  solution  of  the  virtual 
height  integral  for  the  no-field  case  is  known,  we  can  transform  (50)  into  a 
Stilt jes  integral 


(51) 


where  p^ ,  the  no-field  solution,  is  now  used  as  the  integration  variable. 
This  is  a  general  principle  and  can  be  applied  to  a  variety  of  other  models. 
In  the  specific  case  of  a  parabolic  layer  given  by 


2  2  h_hm  2 

fN  =  fc  [1  -  <— 1 '  1 
m 


(52) 


we  have  for  a  half  thickness  of  one  unit  ( y_  =*  1) 

m 


Pn(X) 


Hi 


with  3  =  fc/f 

and  for  the  inverse  function 


(53) 


X  «V4  [2(1  +  32)  -  (1  +  $)2  exp(-2  p  3)  -  ( 1-3)2  exp(2p  3)).  (54) 

n  n 


If  we  divide  pn(1)  -  pft(0)  by  an  even  nutber  ra,  we  obtain  m  intervals  on  the 
Pn  scale.  Now  we  compute  the  corresponding  values  of  X  with  (54),  evaluate 
«*  ST-x  for  those  values  of  x,  and  apply  Simpson's  rule  of  numerical 
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integration  to  obtain  the  parabolic  coefficients. 

Results  from  the  peak  fitting  program  are  shown  in  figure  4.  Again, 
ionograms  were  recorded  in  2-minute  intervals  and  each  set  of  peak  parameters 
shown  was  computed  independently  from  one  single  ionogram. 

ABSORPTION  MEASUREMENTS 

As  mentioned  earlier,  one  data  set  provides,  besides  phases  and  virtual 
height,  eight  amplitudes.  In  most  cases,  the  average  over  eight  amplitudes  is 
significant  in  the  study  of  systematic  variations  (eg,  absorption  as  a 
function  of  frequency  or  at  a  fixed  frequency  as  a  function  of  time).  In 
figure  5,  the  variation  of  the  absorption  with  time  at  a  fixed  frequency  is 
shown  during  a  solar  X-ray  flare  together  with  the  intensity  of  X-ray 
radiation  observed  from  a  satellite.  All  quantities  are  expressed  in  dB,  with 
the  intensity  scale  downwards  for  better  comparison.  It  is  interesting  to 
note  that  at  14h27  the  ionospheric  absorption  is  almost  back  to  normal  for 
this  particular  frequency,  while  the  wideband  and  narrowband  X-ray  intensities 
are  still  above  their  reference  levels  by  20  dB  and  60  dB,  respectively. 
These  values  indicate  a  threshold  which  must  be  exceeded  before  absorption  at 
this  frequency  is  affected.  Approximately  the  same  threshold  values  are 
obtained  if  the  data  are  compared  at  the  onset  of  the  event. 
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CONCLUSIONS 


Some  of  the  algorithms  described  in  this  report  have  been  applied  and 
results  were  published  (Wright,  J.  W.,  1982;  Wright  et  al,  1982),  while  most 
of  the  material  is  new.  It  is  easy  to  see  that  modern  technology  applied  in 
digital  iono sondes  could  open  a  new  area  in  ionospheric  research  and  hf 
communications.  A  few  more  comments  on  the  figures  shown  may  illustrate  this 
better. 

Zenith  distances  of  10°  and  more  (as  shown  in  figure  1)  lasting  for 
approximately  50  minutes  are  expected  to  be  a  regular  feature  during  the  sun¬ 
rise  period.  Similar  effects  must  exist  for  the  sunset  period,  probably 
lasting  longer  with  somewhat  smaller  zenith  angles.  On  the  one  hand,  more 
detailed  systematic  measurements  of  the  apparent  echo  location,  together  with 
the  Doppler  velocity,  could  give  more  insight  into  the  ionization  and 
recombination  process  in  the  F-region.  On  the  other  hand,  the  determination 
of  electron  density  profiles  and  the  height  of  the  F-region  maximum  becomes  an 
entirely  new  and  more  complicated  problem  under  those  conditions. 

Local  F-region  tilts  corresponding  to  zenith  angles  of  approximately  10* 
have  a  significant  effect  on  hf  propagation  and  cannot  be  neglected  in  an  hf 
prediction  model.  Any  attempt  of  hf  direction  finding  or  trananitter  location 
will  fail  under  those  conditions. 
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Sunrise  and  sunset  tilts  could  be  measured  and  become  predictable,  and 


corresponding  corrections  could  be  applied  to  hf  prediction  models.  Tilts  of 
comparable  magnitude  caused  by  gravity  waves  can  appear  at  any  time,  however, 
and  can  have  similar  effects  on  hf  propagation  conditions.  The  periodicities 
of  these  gravity  waves  are  usually  somewhat  shorter  (10-30  minutes;  see 
insert  of  figure  3)  than  the  duration  of  sunrise  and  sunset  tilts  (40-60 
minutes).  Gravity  waves  could  become  predictable,  at  least  in  a  statistical 
sense,  if  correlations  with  other  geophysical  parameters  could  be  established. 
It  should  be  noted  that  even  weak  gravity  waves  can  cause  signal  distortion 
and  fading  in  propagation  through  those  areas,  and  that  the  quality  of  signals 
is  inversely  related  to  gravity  wave  activity. 

The  increased  sensitivity  of  modem  digital  ionosondes  is  best  illus¬ 
trated  in  figure  3.  Changes  in  the  MUF(3000)  on  the  order  of  10  kHz  are 

-4 

already  significant,  corresponding  to  a  relative  accuracy  of  the  order  5:10 
for  an  individual  measurement,  while  the  propagation  factor  M  (3000)  obtained 
from  analog  ionograms  is  typically  accurate  to  1.5*10  for  monthly  median 
values.  This  high  accuracy  may  not  be  important  for  hf  propagation,  but 
defintely  increases  the  detectability  of  natural  and  artifical  disturbances  by 
more  than  one  order  of  magnitude,  provided  that  the  ionogram  sequence  is  fast 
enough  (at  least  one  ionogram  every  5  minutes). 

Similar  comments  can  be  made  about  the  data  presented  in  figure  4.  As 
mentioned  earlier,  the  absolute  value  of  h^  may  be  influenced  by  the  FI-layer, 
if  present.  The  periodicities  and  the  magnitudes  of  the  temporal  variations 
visible,  however,  should  be  accurate.  More  tests  with  simulated  data  may  be 
required. 
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The  length  of  the  frequency  Interval  in  which  ordinary  and  extraordinary 
components  are  overlapping  depends  for  a  given  ionosphere  on  the  pulse  length, 
and  can  easily  be  on  the  order  of  half  the  frequency  range  where  echoes  are 
observed.  The  classical  approach  to  minimize  this  problem  is  to  keep  the 
pulse  length  short.  of  course,  a  shorter  pulse  duration  means  a  wider  spec¬ 
trum,  which  requires  a  wider  bandwidth  of  the  receiver  and  more  power  to  main¬ 
tain  a  reasonably  high  signal- to-noise  ratio.  It  also  causes  poorer  accuracy 
in  parts  of  the  ionogram  where  dispersion  becomes  significant  (Paul,  A.  K. , 
1979) .  Wider  bandwidth  and  higher  power  of  the  transmitted  pulse  causes  more 
interference  to  other  hf  users.  In  general,  a  pulse  duration  of  60-100  l»s 
has  become  standard.  with  the  solution  of  the  o-x  overlap  problem,  the 
pulse  length  could  be  increased,  the  pulse  power  reduced,  and  the  receiver 
bandwidth  made  narrower  without  loss  of  information  at  least  for  F-region 
observations.  For  E-region  observations,  where  often  small  Es  patches  are 
embedded  in  the  regular  E-region,  the  two- wave  solution  (overlap  of  two  echoes 
with  equal  polarization)  would  have  to  be  incorporated  in  the  analysis  program 
if  the  pulse  length  were  increased. 

Besides  possible  improvements  of  the  instrument,  the  separation  of 
ordinary  and  extraordinary  components  as  demonstrated  in  figure  2.  provides  a 
complete  echo  trace  for  each  component  with  high  accuracy,  provided  no  other 
disturbance  is  present  (spread  F,  blanketing  Es)  •  This  is  important  for  al¬ 
most  every  analysis  of  vertical  sounding  data,  like  electron  density  profiles, 
gravity  waves  propagating  through  the  F-region  as  seen  in  the  angle  of 
arrival,  or  virtual  heights  at  fixed  frequencies  as  a  function  of  time, 
ionospheric  drift  measurements,  etc. 
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The  effect  of  increased  absorption  caused  by  solar  flares,  as  shown  in 
figure  5,  is  well  known,  but  few  observations  were  made  measuring  the  fre¬ 
quencies  and  time  dependence  of  the  absorption  and  comparing  it  with  X-ray 
data. 


More  observations  of  this  kind  could  at  least  provide  a  prediction  of  the 
duration  of  such  an  absorption  event  as  a  function  of  frequency,  depending  on 
the  intensity  of  an  observed  flare,  as  long  as  flares  themselves  cannot  be 
predicted. 

One  major  advantage  of  modern  digital  ionosondes  is,  of  course,  the  auto¬ 
mation  of  the  data  analysis.  If  a  minicomputer  is  part  of  the  system,  the 
analysis  could  be  done  in  real  time.  This  can  be  important,  where  real-time 
information  about  the  ionosphere  is  needed  (rocket  experiments  or  timing  of 
artificial  disturbances) .  At  the  same  time,  real-time  processing  could 
eliminate  the  need  of  storing  the  raw  data  by  saving  the  analyzed  results 
only,  which  would  result  in  a  significant  data  compression. 

It  may  be  added  that  the  savings  accomplished  by  the  automatic  analysis 
more  than  compensates  for  the  higher  cost  of  a  modern  digital  system.  In 
addition,  such  ionosondes  can  significantly  improve  our  knowledge  about  the 
ionosphere,  its  physics,  and  its  use  as  a  propagation  median. 
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Figure  1 .  Oblique  echoes  (angle  of  arrival),  Doppler  velocity  and  MUF  (3000) 
deviation  during  sunrise  period. 
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Figure  2.  Separation  of  ordinary  (•)  a.  extraordinary  (x)  echoes  from 
mixed  (M)  echoes. 
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Figure  3.  High-retolution  temporal  variation  of  MUF  (3000). 
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Figure  5 .  Comparison  of  absorption  with  X-ray  flux  during  a  solar  flare. 
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